import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit, least_squares
from scipy.interpolate import splrep, splev

# 设置字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'

def calculate_r_squared(y_actual, y_pred):
    # 计算均值
    y_mean = np.mean(y_actual)

    # 计算残差平方和（RSS）和总平方和（TSS）
    rss = np.sum((y_actual - y_pred) ** 2)
    tss = np.sum((y_actual - y_mean) ** 2)

    # 计算R方
    r_squared = 1 - (rss / tss)

    return r_squared

# 定义拟合函数 V1 = a + b*exp(C1) + c*N1 + d*T
def fit_function(t, a, b, c, d):
    C1, N1, T = t  # t 作为输入的三个变量: C1, N1, T
    return a + b * C1 * np.exp(c * N1) + d * T

# 读取 Excel 文件中的数据
file_path = './assets/normalized_data_no_panimic.xlsx'  # Excel 文件路径
data = pd.read_excel(file_path)

# 提取数据列
C1 = data['C1'].values  # 特征 C1
N1 = data['N1'].values  # 特征 N1
T = data['T'].values  # 特征 T
V1 = data['V1'].values  # 特征 V1

# 定义一个残差函数：实际值与拟合值之差
def residuals(params, t, actual_values):
    return fit_function(t, *params) - actual_values

# 设置初始参数猜测值
initial_params = [1, 1, 1, 1]  # a, b, c, d 的初始值

t = (C1, N1, T)
# 使用 least_squares 进行拟合
result = least_squares(residuals, initial_params, args=(t, V1))

# 提取拟合参数
a, b, c, d = result.x
print(f'拟合参数: a={a:.4f}, b={b:.4f}, c={c:.4f}, d={d:.4f}')

# 计算拟合值
V1_fit = fit_function(t, a, b, c, d)

# 计算R方
r_squared = calculate_r_squared(V1, V1_fit)
print(f'R²值: {r_squared:.4f}')

# 使用样条插值来平滑拟合曲线
spl = splrep(T, V1_fit, s=0)  # s=0 保证插值通过所有数据点
T_smooth = np.linspace(min(T), max(T), 500)  # 平滑的T值，增加采样点数
V1_smooth = splev(T_smooth, spl)  # 平滑后的V1值

# 使用 seaborn 绘图
sns.set_theme(style="darkgrid")
plt.figure(figsize=(10, 6))
sns.scatterplot(x=T, y=V1, color=sns.color_palette("PuBu")[2], label='Actual Data')
sns.lineplot(x=T_smooth, y=V1_smooth, color=sns.color_palette("PuBu")[4], label='Fitted Curve')
plt.ylabel('V(t)', fontsize=14, fontname='Times New Roman')
plt.title('V(t) Model Fit', fontsize=16, fontname='Times New Roman')
plt.legend()
plt.xticks(ticks=np.arange(T.min(), T.max() + 1, 1), fontname='Times New Roman')  # 设置横轴刻度为整数
plt.yticks(fontname='Times New Roman')
plt.show()